# -*- coding: utf-8 -*-
# Scikit-learn-4.py
"""
实验题目4：主成分分析(PCA)降维与可视化
任务描述：
使用主成分分析(PCA)技术对高维葡萄酒数据集进行降维，探索数据的内在结构，并比较降维前后分类模型的性能。

数据集：
使用Scikit-learn内置的葡萄酒数据集（load_wine），包含13个化学特征，用于分类3种葡萄酒。

步骤：

1. 加载葡萄酒数据集，了解数据的基本结构和特征维度
2. 对数据进行标准化处理（使用StandardScaler）
3. 应用PCA进行降维，通过方差解释率确定最佳主成分数量
4. 可视化降维效果：绘制主成分的方差解释率曲线和累积方差解释率曲线
5. 使用前两个主成分绘制数据散点图，观察类别分离情况
6. 比较降维前后使用不同分类器（如SVM、随机森林）的性能差异
7. （可选）使用3D可视化展示前三个主成分的数据分布
"""

# 导入所需库
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import os

# 设置中文字体支持
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("实验四：主成分分析(PCA)降维与可视化")
print("=" * 60)

# 1. 加载葡萄酒数据集
print("1. 加载葡萄酒数据集...")
wine = load_wine()
X = wine.data
y = wine.target
class_names = wine.target_names
feature_names = wine.feature_names

# 特征名称中英文对照
feature_names_en_cn = {
    'alcohol': '酒精 (alcohol)',
    'malic_acid': '苹果酸 (malic_acid)',
    'ash': '灰分 (ash)',
    'alcalinity_of_ash': '灰分碱度 (alcalinity_of_ash)',
    'magnesium': '镁 (magnesium)',
    'total_phenols': '总酚 (total_phenols)',
    'flavanoids': '黄酮类 (flavanoids)',
    'nonflavanoid_phenols': '非黄酮类酚 (nonflavanoid_phenols)',
    'proanthocyanins': '原花青素 (proanthocyanins)',
    'color_intensity': '颜色强度 (color_intensity)',
    'hue': '色调 (hue)',
    'od280/od315_of_diluted_wines': '稀释葡萄酒的OD280/OD315比值 (od280/od315_of_diluted_wines)',
    'proline': '脯氨酸 (proline)'
}

# 创建中英文混合的特征名称列表
feature_names_mixed = [feature_names_en_cn.get(name, name) for name in feature_names]

print(f"数据集形状: {X.shape}")
print(f"目标类别: {class_names}")
print(f"特征数量: {len(feature_names)}")
print("特征名称 (中英文对照):")
for i, name in enumerate(feature_names_mixed):
    print(f"  {i+1:2d}. {name}")

# 创建数据框以便更好地查看数据
df = pd.DataFrame(X, columns=feature_names_mixed)
df['target'] = y
df['class'] = df['target'].apply(lambda x: class_names[x])

print("\n数据集前5行:")
print(df.head())

# 2. 数据标准化
print("\n2. 数据标准化...")
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
print("数据标准化完成!")

# 3. 应用PCA进行降维
print("\n3. 应用PCA进行降维...")
# 首先拟合所有主成分以查看方差解释率
pca_full = PCA()
X_pca_full = pca_full.fit_transform(X_scaled)

# 确定保留多少主成分能解释95%的方差
cumulative_variance = np.cumsum(pca_full.explained_variance_ratio_)
n_components_95 = np.argmax(cumulative_variance >= 0.95) + 1
print(f"保留 {n_components_95} 个主成分可解释 95% 的方差")

# 使用最佳主成分数量重新拟合PCA
pca = PCA(n_components=n_components_95)
X_pca = pca.fit_transform(X_scaled)
print(f"降维后数据形状: {X_pca.shape}")

# 4. 可视化方差解释率
print("\n4. 可视化方差解释率...")
plt.figure(figsize=(12, 5))

# 方差解释率条形图
plt.subplot(1, 2, 1)
components = range(1, len(pca_full.explained_variance_ratio_) + 1)
plt.bar(components, pca_full.explained_variance_ratio_, alpha=0.7, color='skyblue')
plt.xlabel('主成分 (Principal Components)')
plt.ylabel('方差解释率 (Variance Explained Ratio)')
plt.title('各主成分的方差解释率\n(Variance Explained by Each Principal Component)')
plt.xticks(components)
plt.grid(alpha=0.3)

# 累积方差解释率曲线
plt.subplot(1, 2, 2)
plt.plot(components, cumulative_variance, 'ro-', label='累积方差解释率 (Cumulative Variance)')
plt.axhline(y=0.95, color='green', linestyle='--', alpha=0.7, label='95%方差阈值 (95% Variance Threshold)')
plt.axvline(x=n_components_95, color='red', linestyle='--', alpha=0.7, label=f'最佳主成分数: {n_components_95}')
plt.xlabel('主成分数量 (Number of Principal Components)')
plt.ylabel('累积方差解释率 (Cumulative Variance Explained)')
plt.title('累积方差解释率曲线\n(Cumulative Variance Explained Curve)')
plt.legend()
plt.grid(alpha=0.3)

plt.tight_layout()

# 保存图片
save_dir = "Scikit-Learning/img"
os.makedirs(save_dir, exist_ok=True)
variance_path = os.path.join(save_dir, "pca_variance_explained.png")
plt.savefig(variance_path, dpi=300, bbox_inches='tight')
plt.show()
print(f"方差解释率图已保存为 '{variance_path}'")

# 5. 使用前两个主成分绘制数据散点图
print("\n5. 使用前两个主成分绘制数据散点图...")
plt.figure(figsize=(10, 8))
colors = ['red', 'blue', 'green']

for i, class_name in enumerate(class_names):
    plt.scatter(
        X_pca_full[y == i, 0], 
        X_pca_full[y == i, 1], 
        alpha=0.7, 
        c=colors[i], 
        label=class_name,
        s=50
    )

plt.xlabel(f'第一主成分 (PC1) ({pca_full.explained_variance_ratio_[0]*100:.2f}% 方差)')
plt.ylabel(f'第二主成分 (PC2) ({pca_full.explained_variance_ratio_[1]*100:.2f}% 方差)')
plt.title('葡萄酒数据集PCA降维 (前两个主成分)\n(Wine Dataset PCA - First Two Principal Components)')
plt.legend()
plt.grid(alpha=0.3)

# 保存图片
pca_2d_path = os.path.join(save_dir, "pca_2d_scatter.png")
plt.savefig(pca_2d_path, dpi=300, bbox_inches='tight')
plt.show()
print(f"2D散点图已保存为 '{pca_2d_path}'")

# 6. 比较降维前后分类器性能
print("\n6. 比较降维前后分类器性能...")
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
    X_scaled, y, test_size=0.3, random_state=42
)

# 使用原始数据训练模型
svm_original = SVC(random_state=42)
rf_original = RandomForestClassifier(random_state=42)

svm_original_scores = cross_val_score(svm_original, X_train, y_train, cv=5)
rf_original_scores = cross_val_score(rf_original, X_train, y_train, cv=5)

print("原始数据性能 (5折交叉验证):")
print(f"SVM 平均准确率: {svm_original_scores.mean():.4f} (±{svm_original_scores.std():.4f})")
print(f"随机森林 平均准确率: {rf_original_scores.mean():.4f} (±{rf_original_scores.std():.4f})")

# 使用PCA降维后的数据
X_train_pca = pca.transform(X_train)
X_test_pca = pca.transform(X_test)

svm_pca = SVC(random_state=42)
rf_pca = RandomForestClassifier(random_state=42)

svm_pca_scores = cross_val_score(svm_pca, X_train_pca, y_train, cv=5)
rf_pca_scores = cross_val_score(rf_pca, X_train_pca, y_train, cv=5)

print(f"\nPCA降维后数据性能 (保留{n_components_95}个主成分):")
print(f"SVM 平均准确率: {svm_pca_scores.mean():.4f} (±{svm_pca_scores.std():.4f})")
print(f"随机森林 平均准确率: {rf_pca_scores.mean():.4f} (±{rf_pca_scores.std():.4f})")

# 在测试集上评估最终性能
svm_original.fit(X_train, y_train)
svm_pca.fit(X_train_pca, y_train)

rf_original.fit(X_train, y_train)
rf_pca.fit(X_train_pca, y_train)

y_pred_svm_original = svm_original.predict(X_test)
y_pred_svm_pca = svm_pca.predict(X_test_pca)

y_pred_rf_original = rf_original.predict(X_test)
y_pred_rf_pca = rf_pca.predict(X_test_pca)

print(f"\n测试集性能:")
print(f"SVM 原始数据准确率: {accuracy_score(y_test, y_pred_svm_original):.4f}")
print(f"SVM PCA降维后准确率: {accuracy_score(y_test, y_pred_svm_pca):.4f}")
print(f"随机森林 原始数据准确率: {accuracy_score(y_test, y_pred_rf_original):.4f}")
print(f"随机森林 PCA降维后准确率: {accuracy_score(y_test, y_pred_rf_pca):.4f}")

# 7. 使用3D可视化展示前三个主成分的数据分布
print("\n7. 使用3D可视化展示前三个主成分的数据分布...")
if n_components_95 >= 3:
    from mpl_toolkits.mplot3d import Axes3D
    
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')
    
    for i, class_name in enumerate(class_names):
        ax.scatter(
            X_pca_full[y == i, 0], 
            X_pca_full[y == i, 1], 
            X_pca_full[y == i, 2], 
            alpha=0.7, 
            c=colors[i], 
            label=class_name,
            s=50
        )
    
    ax.set_xlabel(f'第一主成分 (PC1) ({pca_full.explained_variance_ratio_[0]*100:.2f}% 方差)')
    ax.set_ylabel(f'第二主成分 (PC2) ({pca_full.explained_variance_ratio_[1]*100:.2f}% 方差)')
    ax.set_zlabel(f'第三主成分 (PC3) ({pca_full.explained_variance_ratio_[2]*100:.2f}% 方差)')
    ax.set_title('葡萄酒数据集PCA降维 (前三个主成分)\n(Wine Dataset PCA - First Three Principal Components)')
    ax.legend()
    
    # 保存图片
    pca_3d_path = os.path.join(save_dir, "pca_3d_scatter.png")
    plt.savefig(pca_3d_path, dpi=300, bbox_inches='tight')
    plt.show()
    print(f"3D散点图已保存为 '{pca_3d_path}'")
else:
    print("主成分数量不足3个，无法进行3D可视化")

# 输出主成分与原始特征的关系
print("\n主成分与原始特征的关系 (前5个主成分):")
components_df = pd.DataFrame(
    pca_full.components_[:5, :],
    columns=feature_names_mixed,
    index=[f'PC{i+1}' for i in range(5)]
)
print(components_df)

print("\nPCA分析完成!")
print("=" * 60)
print("实验四顺利完成!")
print("=" * 60)